#include "stdafx.h"
#include "MarchingCubes.h"

#pragma comment (lib, "comctl32.lib")
#include <CommCtrl.h>


MCVector LinearInterpolation(MCVector v1, MCVector v2, float alpha)
{
	MCVector v;
	if (v1.val!=v2.val)
	{
		v=v1+(v2-v1)*(alpha-v1.val)/(v2.val-v1.val);
	}
	else
		v=v1;
	return v;
}

//Basic marching cubes algorithm, X, Y, Z denote the 3D space.
// cube ///////////////////////////////////////////////////////////////////////
//    v4----- v5
//   /|      /|
//  v7------v6|
//  | |     | |
//  | |v0---|-|v1
//  |/      |/
//  v3------v2
//  the visiting sequence of indices.
//////////////////////////////////////////////////////////////////////////////
TRIANGLE* MarchingCubesAlgorithm(int X, int Y, int Z, MCVector* vertices, float isoValue, int &numberOfTriangles, HWND hwndParent, HINSTANCE g_hinst)
{
	TRIANGLE* triangles = new TRIANGLE[X*Y*Z/3];//array for triangles, it is easy to manage memory by using 1D array to save 3D data.

	numberOfTriangles = 0;//initialize the number of triangles.

	/*int strideYZ = Y*Z;
	int strideZ = Z;*/
	//////////////////////////////////////
	//Progress Bar
	RECT rcClient;  // Client area of parent window.
	HWND hwndPB;    // Handle of progress bar.
	INITCOMMONCONTROLSEX iccInit;
	int iSpacing;   // Controls the size of the progress bar.
	
	// Initialize the common controls.
	iccInit.dwSize = sizeof(iccInit);
	iccInit.dwICC = ICC_PROGRESS_CLASS;
	InitCommonControlsEx(&iccInit);
	
	GetClientRect(hwndParent, &rcClient); 
	iSpacing = (int)(0.2 * (rcClient.right - rcClient.left));

	//Progress Bar
	hwndPB = CreateWindowEx(NULL,
		PROGRESS_CLASS, TEXT("Progress Bar"), 
		WS_CHILD | WS_VISIBLE
		//| WS_OVERLAPPED | WS_CAPTION
		,
		rcClient.left + iSpacing, 
		(int)((rcClient.bottom - rcClient.top)/2 - iSpacing/2),
		300,
		30, 
		hwndParent, (HMENU)1, g_hinst, NULL);
	int t=(X*3)/100;
	SendMessage(hwndPB,PBM_SETRANGE, 0,MAKELPARAM(0,100));
	SendMessage(hwndPB, PBM_SETSTEP, MAKEWPARAM(1, 0), 0);
	////////////////////////////////////////

	for (int i=0; i<X-1; i++)
	{
		for (int j=0; j<Y-1; j++)
		{
			for (int k=0; k<Z-1; k++)
			{
				//map into a cube. The numbering system for vertices needs a deep consideration
				MCVector verts[8];
				

				verts[0] = vertices[k*X*Y+(i+1)*Y+j];
				verts[1] = vertices[k*X*Y+(i+1)*Y+j+1];
				verts[2] = vertices[(k+1)*X*Y+(i+1)*Y+j+1];
				verts[3] = vertices[(k+1)*X*Y+(i+1)*Y+j];
				verts[4] = vertices[k*X*Y+i*Y+j];
				verts[5] = vertices[k*X*Y+i*Y+j+1];
				verts[6] = vertices[(k+1)*X*Y+i*Y+j+1];
				verts[7] = vertices[(k+1)*X*Y+i*Y+j];

				//derive the cube index by comparing with isoValue
				int index = 0;
				for (int n=0;n<8;n++)
				{
					if (verts[n].val<=isoValue)
					{
						index |= (1<<n);
					}
				}

				//check availability from edgeTable
				if (!edgeTable[index])
					continue;

				//derive points of intersection on edges from edgeTable
				MCVector interVerts[12];
				if (edgeTable[index] & 1) interVerts[0] = LinearInterpolation(verts[0], verts[1], isoValue);
				if (edgeTable[index] & 2) interVerts[1] = LinearInterpolation(verts[1], verts[2], isoValue);
				if (edgeTable[index] & 4) interVerts[2] = LinearInterpolation(verts[2], verts[3], isoValue);
				if (edgeTable[index] & 8) interVerts[3] = LinearInterpolation(verts[3], verts[0], isoValue);
				if (edgeTable[index] & 16) interVerts[4] = LinearInterpolation(verts[4], verts[5], isoValue);
				if (edgeTable[index] & 32) interVerts[5] = LinearInterpolation(verts[5], verts[6], isoValue);
				if (edgeTable[index] & 64) interVerts[6] = LinearInterpolation(verts[6], verts[7], isoValue);
				if (edgeTable[index] & 128) interVerts[7] = LinearInterpolation(verts[7], verts[4], isoValue);
				if (edgeTable[index] & 256) interVerts[8] = LinearInterpolation(verts[0], verts[4], isoValue);
				if (edgeTable[index] & 512) interVerts[9] = LinearInterpolation(verts[1], verts[5], isoValue);
				if (edgeTable[index] & 1024) interVerts[10] = LinearInterpolation(verts[2], verts[6], isoValue);
				if (edgeTable[index] & 2048) interVerts[11] = LinearInterpolation(verts[3], verts[7], isoValue);

				//look up triangle table to find collect triangle vertices
				for (int n=0; triTable[index][n]!=-1; n+=3)
				{
					triangles[numberOfTriangles].v[0]=interVerts[triTable[index][n+2]];
					triangles[numberOfTriangles].v[1]=interVerts[triTable[index][n+1]];
					triangles[numberOfTriangles].v[2]=interVerts[triTable[index][n]];
					triangles[numberOfTriangles].normal=((triangles[numberOfTriangles].v[1] - 
						triangles[numberOfTriangles].v[0]).CrossProduct(triangles[numberOfTriangles].v[2] - 
						triangles[numberOfTriangles].v[0])).Normalize();


					numberOfTriangles++;
				}
			}
		}
		if (i%t==0)
		{
			 SendMessage(hwndPB, PBM_STEPIT, 0, 0); 
			//SendMessage(hwndPB, PBM_SETPOS,k*frames/100,0L);
		}

	}
	TRIANGLE* finalTriangles=new TRIANGLE[numberOfTriangles];
	for (int i=0;i<numberOfTriangles;i++)
	{
		finalTriangles[i]=triangles[i];
	}

	delete []triangles;

	SendMessage(hwndPB, PBM_SETPOS, MAKELPARAM(100, 0), 0); 
	DestroyWindow(hwndPB);

	return finalTriangles;
}

///Marching tetrahedron algorithm, X, Y, Z denote the 3D space.
// tetrahedron ///////////////////////////////////////////////////////////////////////
//                      + 0
//                     /|\
//                    / | \
//                   /  |  \
//                  /   |   \
//                 /    |    \
//                /     |     \
//               +-------------+ 1
//              3 \     |     /
//                 \    |    /
//                  \   |   /
//                   \  |  /
//                    \ | /
//                     \|/
//                      +2
//  the visiting sequence of indices.
//  the visiting sequence of indices.
//////////////////////////////////////////////////////////////////////////////
TRIANGLE* MarchingTetrahedronAlgorithm(int X, int Y, int Z, MCVector* vertices, float isoValue, int &numberOfTriangles, HWND hwndParent, HINSTANCE g_hinst)
{
	TRIANGLE* triangles = new TRIANGLE[X*Y*Z/3];//array for triangles, it is easy to manage memory by using 1D array to save 3D data.

	numberOfTriangles = 0;//initialize the number of triangles.

	/*int strideYZ = Y*Z;
	int strideZ = Z;*/

	//////////////////////////////////////
	//Progress Bar
	RECT rcClient;  // Client area of parent window.
	HWND hwndPB;    // Handle of progress bar.
	INITCOMMONCONTROLSEX iccInit;
	int iSpacing;   // Controls the size of the progress bar.
	
	// Initialize the common controls.
	iccInit.dwSize = sizeof(iccInit);
	iccInit.dwICC = ICC_PROGRESS_CLASS;
	InitCommonControlsEx(&iccInit);
	
	GetClientRect(hwndParent, &rcClient); 
	iSpacing = (int)(0.2 * (rcClient.right - rcClient.left));

	//Progress Bar
	hwndPB = CreateWindowEx(NULL,
		PROGRESS_CLASS, TEXT("Progress Bar"), 
		WS_CHILD | WS_VISIBLE
		//| WS_OVERLAPPED | WS_CAPTION
		,
		rcClient.left + iSpacing, 
		(int)((rcClient.bottom - rcClient.top)/2 - iSpacing/2),
		300,
		30, 
		hwndParent, (HMENU)1, g_hinst, NULL);
	int t=(X*3)/100;
	SendMessage(hwndPB,PBM_SETRANGE, 0,MAKELPARAM(0,100));
	SendMessage(hwndPB, PBM_SETSTEP, MAKEWPARAM(1, 0), 0);
	////////////////////////////////////////

	for (int i=0; i<X-1; i++)
	{
		for (int j=0; j<Y-1; j++)
		{
			for (int k=0; k<Z-1; k++)
			{
				//map into a cube. The numbering system for vertices needs a deep consideration
				MCVector verts[8];
				

				verts[0] = vertices[k*X*Y+(i+1)*Y+j];
				verts[1] = vertices[k*X*Y+(i+1)*Y+j+1];
				verts[2] = vertices[(k+1)*X*Y+(i+1)*Y+j+1];
				verts[3] = vertices[(k+1)*X*Y+(i+1)*Y+j];
				verts[4] = vertices[k*X*Y+i*Y+j];
				verts[5] = vertices[k*X*Y+i*Y+j+1];
				verts[6] = vertices[(k+1)*X*Y+i*Y+j+1];
				verts[7] = vertices[(k+1)*X*Y+i*Y+j];

				MarchingTetrahedron(triangles, verts[0], verts[1], verts[2], verts[6],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[0], verts[1], verts[5], verts[6],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[0], verts[5], verts[4], verts[6],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[0], verts[3], verts[7], verts[6],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[0], verts[2], verts[3], verts[6],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[0], verts[4], verts[7], verts[6],numberOfTriangles,isoValue);

				/*MarchingTetrahedron(triangles, verts[0], verts[1], verts[2], verts[3],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[1], verts[2], verts[6], verts[5],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[0], verts[3], verts[7], verts[4],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[4], verts[5], verts[6], verts[7],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[0], verts[1], verts[5], verts[4],numberOfTriangles,isoValue);
				MarchingTetrahedron(triangles, verts[7], verts[6], verts[2], verts[3],numberOfTriangles,isoValue);*/
			}
		}
		
		if (i%t==0)
		{
			 SendMessage(hwndPB, PBM_STEPIT, 0, 0); 
			//SendMessage(hwndPB, PBM_SETPOS,k*frames/100,0L);
		}
	}
	TRIANGLE* finalTriangles=new TRIANGLE[numberOfTriangles];
	for (int i=0;i<numberOfTriangles;i++)
	{
		finalTriangles[i]=triangles[i];
	}

	delete []triangles;
	SendMessage(hwndPB, PBM_SETPOS, MAKELPARAM(100, 0), 0); 
	DestroyWindow(hwndPB);
	return finalTriangles;
}

//search for triangles in one tetrahedron
void MarchingTetrahedron(TRIANGLE* triangles,
						MCVector v1, MCVector v2,
						MCVector v3, MCVector v4,
						int &numberOfTriangles, float isoValue)
{
	int index = 0;
	
	if (v1.val<=isoValue) index |= (1<<0);
	if (v2.val<=isoValue) index |= (1<<1);
	if (v3.val<=isoValue) index |= (1<<2);
	if (v4.val<=isoValue) index |= (1<<3);
	//check availability from edgeTable
	if (!tetrahedronEdgeTable[index])
		return;

	//derive points of intersection on edges from edgeTable
	MCVector interVerts[6];
	if (tetrahedronEdgeTable[index] & 1) interVerts[0] = LinearInterpolation(v1, v2, isoValue);
	if (tetrahedronEdgeTable[index] & 2) interVerts[1] = LinearInterpolation(v2, v3, isoValue);
	if (tetrahedronEdgeTable[index] & 4) interVerts[2] = LinearInterpolation(v3, v4, isoValue);
	if (tetrahedronEdgeTable[index] & 8) interVerts[3] = LinearInterpolation(v4, v1, isoValue);
	if (tetrahedronEdgeTable[index] & 16) interVerts[4] = LinearInterpolation(v2, v4, isoValue);
	if (tetrahedronEdgeTable[index] & 32) interVerts[5] = LinearInterpolation(v1, v3, isoValue);


	//look up triangle table to find collect triangle vertices
	for (int n=0; tetrahedronTriangles[index][n]!=-1; n+=3)
	{
		triangles[numberOfTriangles].v[0]=interVerts[tetrahedronTriangles[index][n+2]];
		triangles[numberOfTriangles].v[1]=interVerts[tetrahedronTriangles[index][n+1]];
		triangles[numberOfTriangles].v[2]=interVerts[tetrahedronTriangles[index][n]];
		triangles[numberOfTriangles].normal=((triangles[numberOfTriangles].v[1] - 
			triangles[numberOfTriangles].v[0]).CrossProduct(triangles[numberOfTriangles].v[2] - 
			triangles[numberOfTriangles].v[0])).Normalize();


		numberOfTriangles++;
	}
	return;
}
///Marching face algorithm, X, Y, Z denote the 3D space.
// face ///////////////////////////////////////////////////////////////////////
//  v0------v1
//  | .    . | 
//  |   ..   |
//  |  .  .  |
//  v3------v2
//  the visiting sequence of indices.
//////////////////////////////////////////////////////////////////////////////
TRIANGLE* MarchingFacesAlgorithm(int X, int Y, int Z, 
	MCVector* vertices, float isoValue, int &numberOfTriangles, HWND hwndParent, HINSTANCE g_hinst)
	{
	TRIANGLE* triangles = new TRIANGLE[X*Y*Z/3];//array for triangles, it is easy to manage memory by using 1D array to save 3D data.

	numberOfTriangles = 0;//initialize the number of triangles.

	/*int strideYZ = Y*Z;
	int strideZ = Z;*/

	//////////////////////////////////////
	//Progress Bar
	RECT rcClient;  // Client area of parent window.
	HWND hwndPB;    // Handle of progress bar.
	INITCOMMONCONTROLSEX iccInit;
	int iSpacing;   // Controls the size of the progress bar.
	
	// Initialize the common controls.
	iccInit.dwSize = sizeof(iccInit);
	iccInit.dwICC = ICC_PROGRESS_CLASS;
	InitCommonControlsEx(&iccInit);
	
	GetClientRect(hwndParent, &rcClient); 
	iSpacing = (int)(0.2 * (rcClient.right - rcClient.left));

	//Progress Bar
	hwndPB = CreateWindowEx(NULL,
		PROGRESS_CLASS, TEXT("Progress Bar"), 
		WS_CHILD | WS_VISIBLE
		//| WS_OVERLAPPED | WS_CAPTION
		,
		rcClient.left + iSpacing, 
		(int)((rcClient.bottom - rcClient.top)/2 - iSpacing/2),
		300,
		30, 
		hwndParent, (HMENU)1, g_hinst, NULL);
	int t=(X*3)/100;
	SendMessage(hwndPB,PBM_SETRANGE, 0,MAKELPARAM(0,100));
	SendMessage(hwndPB, PBM_SETSTEP, MAKEWPARAM(1, 0), 0);
	////////////////////////////////////////

	for (int i=0; i<X-1; i++)
	{
		for (int j=0; j<Y-1; j++)
		{
			for (int k=0; k<Z-1; k++)
			{
				//map into a cube. The numbering system for vertices needs a deep consideration
				MCVector verts[8];
				

				verts[0] = vertices[k*X*Y+(i+1)*Y+j];
				verts[1] = vertices[k*X*Y+(i+1)*Y+j+1];
				verts[2] = vertices[(k+1)*X*Y+(i+1)*Y+j+1];
				verts[3] = vertices[(k+1)*X*Y+(i+1)*Y+j];
				verts[4] = vertices[k*X*Y+i*Y+j];
				verts[5] = vertices[k*X*Y+i*Y+j+1];
				verts[6] = vertices[(k+1)*X*Y+i*Y+j+1];
				verts[7] = vertices[(k+1)*X*Y+i*Y+j];


				MarchingFace(triangles, verts[0], verts[1], verts[2], verts[3],numberOfTriangles,isoValue);
				MarchingFace(triangles, verts[1], verts[2], verts[6], verts[5],numberOfTriangles,isoValue);
				MarchingFace(triangles, verts[0], verts[3], verts[7], verts[4],numberOfTriangles,isoValue);
				MarchingFace(triangles, verts[4], verts[5], verts[6], verts[7],numberOfTriangles,isoValue);
				MarchingFace(triangles, verts[0], verts[1], verts[5], verts[4],numberOfTriangles,isoValue);
				MarchingFace(triangles, verts[7], verts[6], verts[2], verts[3],numberOfTriangles,isoValue);
			}
		}
		
		if (i%t==0)
		{
			 SendMessage(hwndPB, PBM_STEPIT, 0, 0); 
			//SendMessage(hwndPB, PBM_SETPOS,k*frames/100,0L);
		}
	}
	TRIANGLE* finalTriangles=new TRIANGLE[numberOfTriangles];
	for (int i=0;i<numberOfTriangles;i++)
	{
		finalTriangles[i]=triangles[i];
	}

	delete []triangles;
	SendMessage(hwndPB, PBM_SETPOS, MAKELPARAM(100, 0), 0); 
	DestroyWindow(hwndPB);
	return finalTriangles;
}

//search for triangles in one face
void MarchingFace(TRIANGLE* triangles,
						MCVector v1, MCVector v2,
						MCVector v3, MCVector v4,
						int &numberOfTriangles, float isoValue)
{
	int index = 0;
	
	if (v1.val<=isoValue) index |= (1<<0);
	if (v2.val<=isoValue) index |= (1<<1);
	if (v3.val<=isoValue) index |= (1<<2);
	if (v4.val<=isoValue) index |= (1<<3);
	//check availability from edgeTable
	if (!faceEdgeTable[index])
		return;

	//derive points of intersection on edges from edgeTable
	MCVector interVerts[6];
	if (faceEdgeTable[index] & 1) interVerts[0] = LinearInterpolation(v1, v2, isoValue);
	if (faceEdgeTable[index] & 2) interVerts[1] = LinearInterpolation(v2, v3, isoValue);
	if (faceEdgeTable[index] & 4) interVerts[2] = LinearInterpolation(v3, v4, isoValue);
	if (faceEdgeTable[index] & 8) interVerts[3] = LinearInterpolation(v4, v1, isoValue);
	if (faceEdgeTable[index] & 16) interVerts[4] = LinearInterpolation(v2, v4, isoValue);
	if (faceEdgeTable[index] & 32) interVerts[5] = LinearInterpolation(v1, v3, isoValue);


	//look up triangle table to find collect triangle vertices
	for (int n=0; faceTriangles[index][n]!=-1; n+=3)
	{
		triangles[numberOfTriangles].v[0]=interVerts[faceTriangles[index][n+2]];
		triangles[numberOfTriangles].v[1]=interVerts[faceTriangles[index][n+1]];
		triangles[numberOfTriangles].v[2]=interVerts[faceTriangles[index][n]];
		triangles[numberOfTriangles].normal=((triangles[numberOfTriangles].v[1] - 
			triangles[numberOfTriangles].v[0]).CrossProduct(triangles[numberOfTriangles].v[2] - 
			triangles[numberOfTriangles].v[0])).Normalize();


		numberOfTriangles++;
	}
	return;
}